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Motivated by recent experiments with ultra-cold quantum gases in optical lattices we study the 
decay of the staggered moment in the one-dimensional Fermi-Hubbard model starting from a perfect 
Neel state using exact diagonalization and the iTEBD method. This extends previous work in 
which the same problem has been addressed for pure spin Hamiltonians. As a main result, we 
show that the relaxation dynamics of the double occupancy and of the staggered moment are 
different. The former is controlled by the nearest-neighbor tunneling rate while the latter is much 
slower and strongly dependent on the interaction strength, indicating that spin excitations are 
important. This difference in characteristic energy scales for the fast charge dynamics and the 
much slower spin dynamics is also reflected in the real-time evolution of nearest-neighbor density 
and spin correlations. A very interesting time dependence emerges in the von Neumann entropy, 
which at short times increases linearly with a slope proportional to the tunneling matrix element 
while the long-time growth of entanglement is controlled by spin excitations. Our predictions for 
the different relaxation dynamics of the staggered moment and the double occupancy should be 
observable in state-of-t.he art optical lattice experiments. We further compare time averages of the 
double occupancy to both the expectation values in the canonical and diagonal ensemble, which 
quantitatively disagree with each other on finite systems. We relate the question of thermalization 
to the eigenstate thermalization hypothesis. 

PACS numbers: 03.75.-b 71.10.Pm 37.10.Jk 

I. INTRODUCTION 

The non-equilibrium dynamics of order parameters in 
quenches from ordered into disordered phases and vice 
versa has been the topic of many studies, including work 
on Bose-Einstein condensates BEL bosons defined on 
lattice models j3] and systems with antiferromagnetic or¬ 
der B E]. In quantum magnets, the dynamics of the stag¬ 
gered magnetization is a simple yet non-trivial example 
since the Neel state is never an eigenstate of antiferro¬ 
magnetic Heisenberg models. 

In one spatial dimension, since the spontaneous break¬ 
ing of a continuous symmetry is prohibited, starting from 
a state with perfect Neel order, the staggered magnetiza¬ 
tion is expected to decay to zero under the unitary time 
evolution with a SU(2)-symmetric Hamiltonian. This 
problem has been intensely studied for the spin-1/2 XXZ 
chain HHH and one observes a temporal power-law de¬ 
cay of the staggered magnetization to zero for the XX 
case and indications of an exponential decay to zero in 
the interacting case [6j. The quantum quench dynam¬ 
ics starting from the Neel state has attracted additional 
attention since an exact solution for the long-time asymp¬ 
totic behavior could be obtained exploiting the integra- 
bility of the model Bmm- Therefore, the question of 
whether or not the steady state in this quench problem 
can be described by the generalized Gibbs ensemble [15] 
could be addressed with rigor. 

In the context of condensed-matter experiments, the 
decay of Neel order is related to time-resolved spec¬ 
troscopy with Mott insulators in real materials flUffTH] . 

In experiments with ultra-cold quantum gases, it is of¬ 
ten particularly easy to prepare initial real-space prod¬ 


uct states with a high fidelity, which has been used as 
the starting point in several non-equilibrium studies of 
Hubbard- and Heisenberg-type of models [T5IE5] . The 
particular problem of the decay of Neel order has so far 
been addressed in the non-interacting case in one dimen¬ 
sion [24] (where the initial state is an ideal charge den¬ 
sity wave state of one spin component) and for a two- 
dimensional system [25] . Moreover, the decay of a spin 
spiral has been investigated in a two-component Bose gas 
in the strongly-interacting regime, where it can be de¬ 
scribed by the Heisenberg model, in one and two dimen¬ 
sions [26] . The reverse problem, namely the formation of 
antiferromagnetic order in time-dependent protocols is of 
equal relevance since this may provide a path for study¬ 
ing magnetic order in the quantum regime in ultra-cold 
atomic gas experiments EMU, which has been the goal 
of a series of recent experiments [32l - (36] . For other non¬ 
equilibrium experiments with fermions in optical lattices, 
see EMS]- 

In this work, we study the real-time decay of the 
Neel state in the one-dimensional Fermi-Hubbard model, 
which, first, extends previous studies [5] by incorporating 
charge dynamics and second, is motivated by two related 
recent experiments with fermions in one dimension [24] 
and bosons in two dimensions [25]. The Hamiltonian 
reads 

H = -to y~](cJ_|_i -0 .Cj |0 - + h.c.) + U ^2 riifUii (1) 

i i 

where to is the hopping matrix element, U is the onsite 
repulsion, cj a creates a fermion with spin a =f, 4- on site 
i and rii a = c] a Ci^ a . The initial state is given by 

IV’o) = |..., 4-, - - ■) - 


( 2 ) 
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Consequently, we are at half filling. We use the infinite- 
system size time-evolving-block-decimation (iTEBD) al¬ 
gorithm m to compute the time dependence of sev¬ 
eral observables such as the staggered magnetization, the 
double occupancy, nearest-neighbor correlations and the 
von Neumann entropy (we set H = 1). As a main result, 
we demonstrate that the relevant time scales for the re¬ 
laxation of the double occupancy is set by the inverse of 
the hopping matrix element 1/to while for the staggered 
magnetization and nearest-neighbor spin correlations, 
the dynamics is the slower the larger U is. The differ¬ 
ence in the relaxation dynamics can most clearly be dis¬ 
cerned in the strongly-interacting regime U/to > 4. This 
reflects the existence of two characteristic velocities in 
the low-energy, equilibrium physics of strongly interact¬ 
ing one-dimensional systems, namely the spin and charge 
velocity, related to spin-charge separation [41] . Further¬ 
more, there are fingerprints in the time dependence of 
the entanglement entropy. In general, in global quantum 
quenches, one expects a linear increase of S v n(£) ~ t in 
time [42H44] , In our case, we observe a short-time dy¬ 
namics governed by charge excitations where S v pf ~ tot 
while at longer times S v n oc t/U, suggesting that spin 
excitations are relevant for which the energy scale is the 
magnetic exchange constant J = 4fp /U. 


Furthermore, we analyze the dependency of the dou¬ 
ble occupancy on the post-quench values of U/to and we 
investigate whether the steady-state values are thermal 
or not. The latter is a possible scenario for an integrable 
ID model m- We observe that time averages are close 
to the expectation values in the diagonal ensemble [35] , 
while on the system sizes accessible to exact diagonaliza- 
tion, the expectation values in the diagonal and canon¬ 
ical ensemble are clearly different. In this context, we 
also show that the distribution of eigenstate expectation 
values is in general broad, in contrast to systems that 
are expected to thermalize in the framework of the eigen¬ 
state thermalization hypothesis pfoUTTI . The observation 
of broad eigenstate expectation values of observables in 
our model is similar to those of Refs. [T5) TS] made for 
integrable models of interacting spinless fermions. For 
other recent studies of interaction quantum quenches in 
the one-dimensional Fermi-Hubbard model, see I5UH5I] . 
and for studies of the time evolution starting from a per¬ 
fect Neel state in higher dimensions, see munis]. The 
non-equilibrium dynamics starting from this particular 
state yet combined with a sudden expansion into a homo¬ 
geneous empty lattice has been investigated in Ref. HU. 


The plan of this paper is the following. We provide a 
brief overview over the numerical methods and definitions 
in Sec.jTT] Sect ion [ITT] contains our main results, discussing 
the time evolution of observables and von Neumann en¬ 
tropy, steady-state values, thermalization, and the dy¬ 
namics in the strongly interacting regime. We conclude 


with a summary presented in Sec. IV 


II. NUMERICAL METHODS 


In this work we use two wavefunction-based methods, 
exact diagonalization (ED) and iTEBD, to study non¬ 
equilibrium dynamics in the Fermi-Hubbard model. We 
further use a standard density matrix renormalization 
group code (DMRG) to compute ground-state expecta¬ 
tion values m [59]. 


A. iTEBD 


We use Vidal’s iTEBD algorithm for infinite systems to 
calculate the time evolution of the observables of interest 
starting from the perfect Neel state. This method ap¬ 
proximates the true wave-function by a matrix-product 
state ansatz [60] appropriate for the thermodynamic limit 
and is related to time-dependent density matrix renor¬ 
malization group methods [STI (]2 and TEBD for finite 
systems [63]. We use a Trotter-Suzuki break-up of the 
time-evolution operator with a time step that is cho¬ 
sen small enough to resolve high-frequency oscillations at 
large U/to- The maximum number of states is bounded 
by X max — 1024. We compaied luns with different Xmax 
and show only data for which the results are indistin¬ 
guishable on the scale of the figures. 

Compared to its siblings - the time-dependent den¬ 
sity matrix renormalization group method EUE2] and 
TEBD [63] - the advantage of iTEBD is clearly that it 
is set up directly for the thermodynamic limit. More¬ 
over, both TEBD and iTEBD are particularly well-suited 
for problems in which the initial state has an exact 
matrix-product state representation, which applies to 
our situation. All theses approaches rely on approxi¬ 
mating the time-evolved wave-function through matrix- 
product states which only gives a faithful representation 
if the time-evolved wave-function does not encode a large 
amount of entanglement [60]. While our initial state is 
not entangled, the entanglement in a global quantum 
quench like ours grows linearly in time (see, e.g., [42] ). 
which results in an exponential increase of computational 
effort [60]. Thus, as the time evolution progresses, even¬ 
tually, going from time step t to t + At will consume 
more computational time than the whole previous cal¬ 
culation. By keeping the discarded weight constant in 
every step, one accounts for the time-dependent increase 
of the entanglement entropy, and by carrying out simu¬ 
lations with a different discarded weight the accuracy of 
the data can be controlled. While the linear increase of 
the entanglement entropy with time is generic to a global 
quantum quench, the actual quench, model parameters 
and the observable determine the actual numerical costs 
such that no general prediction of numerical effort and 
accuracy is possible. 
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B. Exact diagonalization 

Our second method is exact diagonalization. We per¬ 
form the time evolution in a truncated Krylov space (see 
[B3t f° r a review and references). To be able to treat 
larger systems we exploit symmetries of the Hamiltonian 
([I]), namely conservation of total particle number N, to¬ 
tal spin S z , invariance under lattice translations (quasi¬ 
momentum k), the parity and spin-flip symmetry. In ED 
simulations, we use periodic boundary conditions, the 
number of sites is denoted by L. 

C. Observables 

Key quantities in our analysis are the double occu¬ 
pancy 

1 L 

= t n 4>> (3) 

i =1 

where the associated operator is d = n it n ii- The 

staggered magnetization is 

1 L 

™s(t) =(n n - n^). (4) 

i=1 

We further study nearest-neighbor density and spin cor¬ 
relations defined as W = (rqrq+i) and 5) = (S/S~ +1 ), 
with rii = nq + 714 and Sf = (nq- — nq.)/2. The von 
Neumann entropy for a central cut through the system 
is computed from 

Svn =-tr [pa In Pa], (5) 

where pa is the reduced density matrix of one half of the 
system. 

III. RESULTS 

A. Time evolution and characteristic time scales 

1. Double occupancy and staggered moment 

Figures [ija) and (b) show the time evolution of the 
double occupancy d(t) and of the staggered magnetiza¬ 
tion m s (t ), respectively, obtained from iTEBD simula¬ 
tions. While the double occupancy rapidly approaches 
a time-independent regime for all values of U/to consid¬ 
ered here, the relaxation of the staggered magnetization 
towards m s = 0 is much slower. It is very instructive to 
replot m s {t) versus t/U [inset in Fig. Qb)]. This results 
in a collapse of the data for U > 4t 0) which is the better 
the larger U/to is. Therefore, the relaxation of double 
occupancy and staggered magnetization occur at differ¬ 
ent time scales 1/to and U, respectively. This suggests 



FIG. 1. (Color online) (a) Double occupancy d(t) and (b) 
staggered magnetization m a (t) as a function of time during 
the quench from the Neel state to U/to = 0,4,8,16 (iTEBD 
data). Dashed lines in (a): expectation value ddiag in the 
diagonal ensemble Eq. < |10| ) from ED (L = 10). Inset in (b): 
m s (t) plotted versus t/U for U/to = 4, 8,16. 

that the relaxation of spin-related quantities is set by the 
magnetic exchange matrix element given by J = 4 t^/U 
for large U/to■ Both quantities further exhibit coherent 
oscillations that decay during the approach to a station¬ 
ary value. For the double occupancy the frequency is 
given by uj = U for large U t 0 . By contrast, the pe¬ 
riod of oscillations in m s (t) increases in the large U/to 
limit. This is expected, since in the Heisenberg limit the 
period of oscillations is 1/(2J) with J = 4C7/[B|. Note 
that the non-interacting case has recently been studied 
comprehensively in |6i5| and that our iTEBD results agree 
with the analytical solution for the U = 0 case [6] [66] [67]. 

The short-time dynamics of both quantities (O repre¬ 
senting an observable) can be obtained analytically by 
expanding the time-evolution operator: 

(0(t)) « {i/o\0\ipo) +i{'/’o\[H 1 6]\i/j 0 )t 

-^o\[H,[H,d}Uo)t 2 + 0(t 3 ). (6) 

For both double occupancy and staggered magnetization, 
the leading time dependence is ~ t 2 and comes from 
(ijjo\HOH\ipo) °c w hich is independent of U. Hence, 
the nontrivial U -dependence cannot be deduced from this 
short-time dynamics. Second-order time-dependent per¬ 
turbation theory in t 0 /U gives 

d (t)=u2 sm (t)’ (7) 

. . 1 8tn . 2 fUt\ /0 . 

’"*«) = 2 -Jp Sm ( t ) < 8 > 
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FIG. 2. (Color online) Time dependence of the staggered 
moment for Heisenberg (solid line) and Hubbard model with 
a large U/to = 16,32 (iTEBD data), plotted versus time 
measured in inverse units of the magnetic exchange constant 
J = 4to/(7. The arrow indicates the small-amplitude oscilla¬ 
tions in the short-time dynamics for finite U /to < oo whose 
frequency is given by U. 


These expressions agree with our numerical data for 
U/t 0 > 16. 



FIG. 3. (Color online) (a) Nearest-neighbor charge correla¬ 
tions Ni(t) = (mm+i) and (b) nearest neighbor, longitudinal 
spin correlations Si(t) = (S/S/ +1 ) as a function of time for 
U/t 0 = 0,4,8,16 (iTEBD results). 


2. Comparison to Heisenberg model 

For completeness, we show that the time dependence 
of the staggered magnetization in the large U/to limit 
approaches the one of the spin-1/2 Heisenberg model 
H = JJ2i S i ■ Sj+1, where J is the magnetic exchange 
coupling. We expect the time evolution of m s (t ) to be 
identical in both models in the limit of large U/J since 
the Heisenberg model is derived from the Fermi-Hubbard 
model via a Schrieffer-Wolff transformation that projects 
onto the subspace of vanishing double occupancy [SB] ■ 

The comparison is shown in Fig. [2j where we present 
iTEBD results for U/to = 16,32 and the pure spin sys¬ 
tem. We see that the results for the two models be¬ 
come quantitatively similar for large U/to. Moreover, the 
short-time dynamics in m s (t ), namely the small initial os¬ 
cillations (see the arrow in the figure), disappear as U/to 
increases, indicating the complete suppression of short- 
time charge dynamics. This is accompanied by a shrink¬ 
ing of the time window in which the short-time dynamics 
is governed by 5m s (t) = m s (t) — m s (t = 0) oc (tot) 2 [see 
Eq. ([6])] which gets replaced by 8m s (t) oc ( Jt) 2 (the latter 
follows from considering the Heisenberg model). 


3. Nearest-neighbor correlations 

The time dependence of nearest-neighbor density cor¬ 
relations Ni(t) [Fig. [3](a)] and spin correlations S*(i) 
[Fig. [3](b)], bears similarities with the one of the double 
occupancy and the staggered moment, respectively. The 
density correlator undergoes a rapid decrease towards a 
stationary state that happens during the first tunneling 
time and then exhibits oscillations with a [/-dependent 
frequency. On the contrary, the relaxation dynamics of 


the spin correlator is much slower, and again controlled 
by U (the data for Si(t) can be collapsed in the U/to > 4 
regime by plotting them versus to/U, analogous to the 
staggered moment). 


4. Von Neumann entropy 


The existence of two different time scales for the re¬ 
laxation dynamics of the double occupancy and the stag¬ 
gered magnetization translates into an interesting time 
dependence of the von Neumann entropy (see Fig. [ 4 ]). At 
short times t < 0.5/to, ~ t with a prefactor that 

is independent of U, while for t > 0.5/to, the time de¬ 
pendence crosses over to a linear increase with a strongly 
[/-dependent slope. Plotting S v n versus t/U results in a 
collapse of the data [see the inset in Fig. [ 4 ], comparable 
to the behavior of the staggered magnetization. 

The prefactor c s of the linear increase of the von Neu¬ 
mann entropy is related to the existence of gapless modes 
and given by the characteristic velocities [44]. We have 
extracted the prefactor of S v n from the increase in the 
[/-dependent regime, shown in Fig. [5] It turns out to 
be a monotonically decreasing function of U/t 0 . We fur¬ 
ther compare c s to the exact value of the spinon velocity 
vf A known from the Bethe ansatz [50K7TI (dashed line 
in Fig. [ 5 ]): 


v 


BA 

s 


hQnto/U) 

°I 0 (27rfo/[/) 


(9) 


(Io and I\ are modified Bessel functions of the first kind). 
Both c s and clearly have a very similar dependence 
on U/to, unambiguously showing that the long-time dy¬ 
namics of the entanglement entropy are controlled by 












5 



FIG. 4. (Color online) Von Neumann entropy S v n for a 
central cut through the system as a function of time for 
U/to = 0,4,8,16 (iTEBD results). Inset: Svn plotted ver¬ 
sus t/U for U/t 0 = 4, 8,16. 

spin excitations. 

B. Time averages of double occupancy 

In the analysis of time averages, it is instructive to 
compare them to the expectation values in the diagonal 
and canonical ensemble. The diagonal ensemble is de¬ 
fined as jT5] 

Odiag = M 2 (<a|0|a), (10) 

Ot 

where |ct) are post-quench eigenstates (H\a) = E a \a)) 
and c a = (V’o|ct) are the overlaps between the initial state 
and post-quench eigenstates. Odiag is the long-time av¬ 
erage of (O) @5], where degeneracies do not enter. 

Given that the double occupancy can routinely be mea¬ 
sured in quantum gas experiments [2U |37j , we concen¬ 
trate the following discussion on this quantity. The val¬ 
ues for ddiag computed for L = 10 using ED are included 
in Fig.JlJa) as dashed lines. Clearly, the time-dependent 
iTEBD data are very close to ddiag and seem to approach 
this value as the amplitude of oscillations decays. 



FIG. 5. (Color online) Characteristic velocities c s extracted 
from the time dependence of the von Neumann entropy SVn 
in the [/-dependent regime t > 0.5/to, plotted versus U/to 
(circles). For comparison we include the exact values vf A 
(dashed line) of the spin velocity known from the Bethe-ansatz 
solution m- 



FIG. 6. (Color online) Finite-size scaling of the expectation 
value of the double occupancy in the diagonal ensemble (cir¬ 
cles) ddiag for (a) U = 4to and (b) U = 8to (circles: ED data 
for L = 4,6,8,10, star: time average d from iTEBD). The 
figure also includes the expectation values in the canonical 
ensemble d can (squares). 

To get a feeling for the system-size dependence, we 
show ddiag versus 1 /L for (a) U = 4t 0 and (b) U = 8 1 0 , 
together with d extracted from iTEBD simulations plot¬ 
ted at 1/L = 0 in Figs. |6ja) and (b), respectively. The 
finite-size dependence of the data for ddiag is consistent 
with d diag (L) —► d as system size increases. We should 
stress, though, that the time average of the double occu¬ 
pancy itself could change if we were able to reach longer 
times with iTEBD. 

The expectation value in the canonical ensemble is 
computed from 

Ocan = tr \pO \, (11) 

where p = exp(—/3H)/Z with Z the partition function, 
all evaluated at fixed N = L and vanishing total spin 
= 0- The temperature T = 1//3 is fixed by 
requiring that 

E = {ipo\H\ipo) = ti[pH\. (12) 

While in our problem E = 0, independently of the post- 
quench value of U, the canonical temperature T clearly is 
a function of U/ J since the post-quench ground-state en¬ 
ergy E gs (U), defining for each U/J the zero-temperature 
reference point, depends on U/J. To illustrate this point, 
we introduce the excess energy 

SE = E-E gs (U). (13) 

The canonical temperature T/U expressed in units of 
U and the excess energy SE are plotted versus U/J in 
the main panel and inset of Fig. [TJ respectively. Both 
TjU and SE are monotonously increasing functions of 
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FIG. 7. (Color online) Canonical temperature T/U (main 
panel) and excess energy SE (see Eqs. (121 and ©) versus 
U/J for L w 10 (ED results). 



FIG. 8. (Color online) (a) Time averages (circles) of the dou¬ 
ble occupancy as a function of U/to- Time averages are ob¬ 
tained by averaging over full periods of the oscillations. Tri¬ 
angles denote the ground-state expectation values (computed 
with DMRG for L = 64 and open boundary conditions) for 
comparison, stars are the expectation values in the canoni¬ 
cal ensemble Eq. © computed with ED for L = 10. (b) 
Relative difference between canonical and diagonal ensemble 

Ad re l — (ddiag dean )/ddiag- 


as anticipated from Fig. [lja), d ~ ddi ag for the accessi¬ 
ble time scales/system sizes [data for ddiag not shown in 
Fig- Ufa)]. 

Second, the time averages are above the ground- 
state expectation values. This behavior is, in the large 
U/to limit, somewhat unexpected, given the known non¬ 
monotonic temperature dependence of d. As a function of 
T, the equilibrium double occupancy d{T) first decreases 
from its zero-temperature value and then increases for 
large T towards d(T = oo) = 1/4 (see {73j [74]). The 
position of the minimum in d(T ) can be interpreted as a 
scale for the separation of the spin- versus charge exci¬ 
tation dominated temperature regime. Since we do not 
observe d < d gs up to U/to = 64, we conclude that the 
initial state always mixes in doublons from the upper 
Hubbard band and not just the virtual doublons present 
in the ground state. For the accessible system sizes, this 
is confirmed by the discussion presented in Sec. |HICl 


We further observe the known d gs oc 1/C/ 2 behavior [751 
[75] in the large U/to regime (also obeyed by d). The value 
of d = 1/4, which is the infinite-temperature expectation 
value at U = 0, is approached by d and ddiag as U/to is 
lowered (see Fig. [8]). 

Since the system is integrable, it is not surprising that 
the expectation values in the canonical ensemble are dif¬ 
ferent from the ones in the diagonal ensemble. The 
canonical ensemble has been computed for a small sys¬ 
tem using exact diagonalization, and therefore, a quan¬ 
titative comparison only makes sense by comparing to 
the diagonal ensemble but not to the iTEBD time av¬ 
erages. The relative difference is shown in Fig. [s][b) for 
L = 10 and can be quite large. At least for the accessible 
system sizes (see Fig. 6]), this difference does not seem 
to become smaller. Therefore, we do not observe ther- 
malization in this model for the quench protocol studied 
here. Nonetheless, the qualitative dependence of d , d can 
and ddi a g on U/to is quite similar. 


C. Connection to eigenstate thermalization 
hypothesis 


U/J as U/J is lowered. At U/J = oo, SE is zero since 
the initial state is in the ground state manifold in that 
limit. As U/J decreases, E = 0 moves towards the mid¬ 
dle of the many-body spectrum (see also the discussion in 
Sec. IIIC2 and Fig. |9j) and eventually, at U = 0, it trans¬ 
lates into an infinite temperature (see also [72]). It is thus 
more appropriate to express T in units of U rather than 
J in the large U/J regime since this results in T/U —l 0 
for U/J ^ oo [71. 

From the time-dependent data shown in Fig. 0a), we 
extract the time averages d of the double occupancy. 
These are displayed in Fig. |8][a) versus U/to (circles) to¬ 
gether with the expectation values d gs in the ground state 
(triangles, DMRG data) and the expectation value d can 
in the canonical ensemble (stars). First, we observe that, 


1. Eigenstate expectation values 

One popular framework to understand thermalization 
in closed many-body systems is the eigenstate thermal¬ 
ization hypothesis (ETH) [451447] . It states that Odiag = 
O mc , where O mc is the expectation value in the micro- 
canonical ensemble, if the expectation values (a\0\a) of 
O (a local observable) in post-quench eigenstates only de¬ 
pend on energy E in the thermodynamic limit (the latter 
also assuming a narrow initial state 03172]). In other 
words, expectation values computed in a typical many- 
body eigenstate (which should be the vast majority of all 
states) already yield thermal behavior. For sufficiently 
large systems, expectation values in the microcanonical 
and canonical ensemble should agree with each other. 

On a finite system accessible to exact diagonalization, 
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FIG. 9. (Color online) Post-quench eigenstate expecta¬ 
tion values for the double occupancy for various interaction 
strengths: (a) U/to = 0, (b) U/to = 4, (c) U/to = 8, (d) 
U/to = 16 (ED data for L = 10). The vertical dashed lines 
mark the quench energy E = 0 for our initial state. The inset 
in (d) shows a blow-up of the first doublon band. The Neel 
state is doubly degenerate (denoted by |-0o) an d |V’o}) and the 
linear combinations |i/j±) = (l^o) ± \‘4>o))/V^ live in the total 
quasi-momentum k = 0 and k = n subspaces. 


validity of the ETH manifests itself in a narrow width 
of {a\0\a) at a fixed energy E for a generic quantum 
system, while for a ID integrable system, (a|0|a) can 
be very broad for a given energy, due to the existence of 
many non-trivial (local) conservation laws resulting in a 
large fraction of degeneracies. The picture has been stud¬ 


ied and often verified (see, e.g., H2 M HO EH C3EB), 
the important question being how quickly the distribu¬ 
tions of (a|0|a) become sufficiently narrow as system size 
increases. Recent work suggests that for a generic sys¬ 
tem, this is exponentially fast in L ([52], see also C2!S3), 
while for an integrable system, the decay of the width of 
{a\0\a) at a given E is at most power-law .52] [MJ [55] 
(see also [55]). 

Here, we exclusively analyze the distribution of post- 
quench eigenstate-expectation values of the double occu¬ 
pancy. These are presented in Figs. |9]]a)-(d) for U/to = 
0,4,8,16. For U/to > 4, the distributions have a very 
regular structure inherited from the U/to = oo limit, 
where the double occupancy is a conserved quantity. 
There is one band for each possible value of (a|d|a) 
(for the parameters of the figure, L = 10 these are 
L(a\d\a) = 0,1, 2,3,4, 5). For a nonzero and small t 0 /U, 
the exact degeneracy in these bands is lifted while the 
structure as such is preserved on these small systems. In 
the lowest band, the effect of to 7 ^ 0 is to lower the energy 
from the degenerate U/to = oo ground-state manifold 
at E = 0 towards the correlated ground state, result¬ 
ing at the same time in an increase of (a|(i|a) towards 
its nonzero ground-state expectation value. This lowest 
band is very sharp and its negative slope translates into 
the decrease of d = d(T) from its zero-temperature value 
as a function of temperature at low T mm, which per¬ 
sists as long as the dL = 0 band remains well separated 
from the dL = 1 band. 

At smaller U/to, the bands eventually start to overlap 
and they become very broad at a fixed energy (compare 
the discussion in HU [57] for other models). At U = 0, 
the distribution of (a|d|a) becomes flat, resulting in an 

essentially energy-independent mean value of (a|d|a) ~ 
1/4. 


2. Properties of the specific initial state 


Our initial state has a mean energy of E = 0 
and a width (in the diagonal ensemble) of Odiag = 
-y/ (ipo\H 2 \4’o) = toALL, which is independent of U. This 
is indicated by the shaded areas in Fig. [9] For large 
U/to , primarily the very narrow first band is sampled and 
E = 0 sits at the high-energy edge of the first, dL = 0 
band (recall that for U/to = °°, dL takes integer val¬ 
ues). Therefore, the initial state asymmetrically mixes 
in eigenstates with too large values of (a|d|a) from first, 
states in the dL = 0 band at E < 0 and second, from the 
band with dL = 1 (the latter follows from analyzing the 
distribution of |c a | 2 ). Hence, the overall structure of the 


distribution of (a|d|a) combined with the distribution of 
| Cq, | 2 is consistent with t he obs ervation that d > d can at 
large U/to (compare Sec. Ill B). 

At very small U/to, the initial state samples the bulk 
of the system where the density of states is large. At 
U = 0, the corresponding canonical temperature derived 

















from the quench energy is infinite and since (a|d|a) does 
not depend much on energy, we must find d = ddiag = 
4an —* 1/4 as L increases, consistent with the discussion 
in Sec. IIIB At intermediate U/to, the initial state sam¬ 
ples several overlapping and partially very broad bands 
of the (a\d\a) distribution [see, e.g., the case of U/to = 4 
shown in Fig. [9j[b)] . Therefore, based on the structure 
of the eigenstate-expectation-value distributions at the 
quench energy, we expect deviations between thermal be¬ 
havior at intermediate and large U / J, consistent with 
our previous analysis. In conclusion, we stress that the 
quench energy alone is not a sufficient criterion for the 
analysis of finite-system size data, but that the actual 
distribution of overlaps |cq,| 2 crucially determines which 
bands are involved (see also the discussion in 72]). 


IV. SUMMARY AND CONCLUSION 


In this work, we studied the relaxation dynamics in the 
one-dimensional Fermi-Hubbard model starting from a 
perfect Neel state as a function of the interaction strength 
U/to■ As a main result, we reported evidence that the 


relaxation dynamics of the staggered moment, spin cor¬ 
relations and of the von Neumann entropy at long times 
is controlled by spin excitations, while the double oc¬ 
cupancy undergoes a much faster dynamics controlled 
by charge excitations. The slope c s of the increase of 
the von Neumann entropy S v n = c s t is very similar to 
the exact spinon velocity known from the Bethe ansatz. 
This separation of time scales for double occupancy ver¬ 
sus staggered magnetization could be accessible in state- 
of-the-art quantum gas experiments. 

We further demonstrated that the time averages of 
the double occupancy are different from the expecta¬ 
tion values in the canonical ensemble. Nonetheless, both 
quantities exhibit the same qualitative dependence on 
U/to- Finally, we made a connection to the eigenstate 
thermalization hypothesis by showing that the eigen¬ 
state expectation values of the double occupancy are, 
in general, broadly distributed with no well-defined de¬ 
pendence on energy only, characteristic for an integrable 
one-dimensional system. 
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